#include "MethodMinimization.h"

using namespace FitMethods;
using namespace RooFit;
using namespace Utility;

FitResult MethodMinimization::runAnalysisCase(AnalysisConfig config, FitInput input) {
    if(input.getSignalTemplates().size()!=1 || input.getBackgroundTemplates().size()<1) {
        throw FitMethodException("MethodMinimization needs exactly one signal template (because of subtrBgdFromData) and at least one background template, but was called with "+boost::lexical_cast<string>(input.getSignalTemplates().size())+" signal templates and "+boost::lexical_cast<string>(input.getBackgroundTemplates().size())+"background templates");
    }
    RooRealVar * isol = new RooRealVar("isol","E_{T}^{cone40}/E_{T}",-0.2,1.5);
    gc->add(isol);
    RooRealVar & nSignal = *new RooRealVar("nSignal","amount of signal events", 500, 0, 1000000000);
    gc->add(&nSignal);
    vector <RooRealVar *> signalVector;
    signalVector.push_back(&nSignal);
    vector <RooRealVar *> backgroundVector;
    for(int i = 0;i<(int)input.getBackgroundTemplates().size();i++) {
        RooRealVar & bgdTempl = *new RooRealVar((string("nBackground")+boost::lexical_cast<string>(i)).c_str(),"amount of signal events", 500, 0, 1000000000);
        gc->add(&bgdTempl);
        backgroundVector.push_back(&bgdTempl);
        //WARNING: irregular behavior
        if(i==0) {
            Utility::subtrBgdFromData(input.getData(),input.getBackgroundTemplates()[i],*backgroundVector[i],nSignal,config.getIsolThr());
        } else {
            RooRealVar temporaryVariable("nTemporary","used for subtrBgdFromData", 0, 0, 1000000000);
            Utility::subtrBgdFromData(input.getData(),input.getBackgroundTemplates()[i],*backgroundVector[i],temporaryVariable,config.getIsolThr());
        }
        cout << "Starting nbgd[" << i << "]: " << backgroundVector[i]->getVal() << endl;
        cout << "Starting nSignal: " << nSignal.getVal() << endl;
    }
    for(int i = 0;i<(int)backgroundVector.size();i++) if(backgroundVector[i]->getVal()<3) backgroundVector[i]->setVal(3);
    for(int i = 0;i<(int)signalVector.size();i++) if(signalVector[i]->getVal()<3) signalVector[i]->setVal(3);
    vector <RooDataHist *> bgdDataHistVector = createRooDataHistVector("background",input.getBackgroundTemplates(),*isol,gc);
    vector <RooHistPdf *> bgdHistPdfVector = createRooHistPdfVector("background",bgdDataHistVector,*isol,gc);
    vector <RooDataHist *> sigDataHistVector = createRooDataHistVector("signal",input.getSignalTemplates(),*isol,gc);
    vector <RooHistPdf *> sigHistPdfVector = createRooHistPdfVector("signal",sigDataHistVector,*isol,gc);
    RooAddPdf * model = createRooPdfAsSum("model",sigHistPdfVector,bgdHistPdfVector,signalVector,backgroundVector,gc);
    /// WARNING: Isolation used improperly
    RooDataHist * rooDataHist = new RooDataHist("rooDataHist", "Isolation", *isol, input.getData());
    gc->add(rooDataHist);
    string name;
    if(config.getIsMethodEXLL()) {
        name="Extended Likelihood Method";
        model->fitTo(*rooDataHist,Extended(kTRUE),SumW2Error(kFALSE));
    }
    else {
        name="Chi2 Minimization Method";
        model->chi2FitTo(*rooDataHist,Extended(kTRUE),SumW2Error(kFALSE));
        /*
        RooChi2Var chi2("chi2var","chi2var",*model,*rooDataHist,Extended(kTRUE),SumW2Error(kFALSE));
        RooMinuit minimizer(chi2);
        minimizer.migrad();
        minimizer.hesse();
        minimizer.minos();
        */
    }
    Utility::correctToRange(nSignal, input.getSignalTemplates()[0], config.getIsolThr(), true);
    FitResult res(name,signalVector,backgroundVector); //< The result object, to be filled and returned

    /** Chi square / ndof goodness of fit calculation */
    RooPlot* frame = isol->frame();
    rooDataHist->plotOn(frame);
    model->plotOn(frame, LineColor(kBlue));
    res.setChiSquare(frame->chiSquare(input.getSignalTemplates().size()+input.getBackgroundTemplates().size()));
   
    /** Setting plausible values for integral below isolation threshold */
    vector <RooRealVar *> nBgdIsolThVector, nSigIsolThVector;
    /// WARNING: multiple background vectors ignored for nBgdIsolTh
    nBgdIsolThVector.push_back(backgroundVector[0]);
    nSigIsolThVector.push_back(signalVector[0]);
    res.setBgdBelowIsolTh(nBgdIsolThVector);
    res.setSigBelowIsolTh(nSigIsolThVector);
    return res;
}
